Last updated: 2024-08-20

Checks: 7 0

Knit directory: DCM_snRNAseq/

This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20240606) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version c979c9e. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/figure/

Untracked files:
    Untracked:  analysis/omnipathr-log/
    Untracked:  code/Add_metadata.R
    Untracked:  code/Check DYSF.R
    Untracked:  code/Clustering_genes.R
    Untracked:  code/DE_5_percent.Rmd
    Untracked:  code/DE_5_percent.html
    Untracked:  code/DE_5_percent/
    Untracked:  code/DE_CM_test1.R
    Untracked:  code/DE_no_401.Rmd
    Untracked:  code/DE_no_401.html
    Untracked:  code/DE_no_401/
    Untracked:  code/Differential abundance test1.R
    Untracked:  code/Differential_Expression_edgeR_All.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_2.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_2.html
    Untracked:  code/Differential_Expression_edgeR_All_Age.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_Age_2.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_Age_2.html
    Untracked:  code/Differential_Expression_edgeR_All_groups.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_groups.html
    Untracked:  code/Differential_Expression_edgeR_All_groups_2.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_groups_2.html
    Untracked:  code/Differential_Expression_edgeR_All_groups_3.Rmd
    Untracked:  code/Differential_Expression_edgeR_All_groups_3.html
    Untracked:  code/PCA and CCA analysis test.R
    Untracked:  code/Pseudobulk DCM HC.R
    Untracked:  code/QC_integration_annotation_orig.Rmd
    Untracked:  code/UpSet_plot_DEGs.Rmd
    Untracked:  code/Volcano_highlighted_genes.R
    Untracked:  code/ezInteractiveTable.Rmd
    Untracked:  code/ezInteractiveTable.html
    Untracked:  code/old.R
    Untracked:  core
    Untracked:  data/Cellbender_output/
    Untracked:  data/Cellranger_output/
    Untracked:  data/DCM_Clinical_data.xlsx
    Untracked:  data/DCM_Clinical_data_26.xlsx
    Untracked:  data/Raw/
    Untracked:  omnipathr-log/
    Untracked:  output/Cardiomyocytes_DA_DE/
    Untracked:  output/Cardiomyocytes_DA_DE_26/
    Untracked:  output/Cardiomyocytes_subclustering/
    Untracked:  output/Cardiomyocytes_subclustering_26/
    Untracked:  output/Differential_expression_edgeR_All/
    Untracked:  output/Differential_expression_edgeR_All_Age/
    Untracked:  output/QC_integration_annotation/
    Untracked:  output/QC_integration_annotation_26/

Unstaged changes:
    Modified:   analysis/Cardiomyocytes_DA_DE.Rmd
    Deleted:    analysis/Differential_Expression_edgeR_All.Rmd
    Deleted:    analysis/Differential_Expression_edgeR_All_Age.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the repository in which changes were made to the R Markdown (analysis/Cardiomyocytes_DA_DE_26.Rmd) and HTML (docs/Cardiomyocytes_DA_DE_26.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File Version Author Date Message
Rmd c979c9e GinoBonazza 2024-08-20 wflow_publish("analysis/Cardiomyocytes_DA_DE_26.Rmd")
html 182f9de GinoBonazza 2024-08-19 Build site.
html 5cd0f74 GinoBonazza 2024-08-01 Build site.
Rmd f3f2075 GinoBonazza 2024-08-01 wflow_publish("analysis/Cardiomyocytes_DA_DE_26.Rmd")
html 61e8d17 GinoBonazza 2024-08-01 Build site.
Rmd 9e4b482 GinoBonazza 2024-08-01 wflow_publish("analysis/Cardiomyocytes_DA_DE_26.Rmd")

Setup

# Get current file name to make folder
current_file <- "Cardiomyocytes_DA_DE_26"

# Load libraries
library(here)
library(readr)
library(readxl)
library(xlsx)
library(Seurat)
library(DropletUtils)
library(Matrix)
library(scDblFinder)
library(scCustomize)
library(dplyr)
library(ggplot2)
library(magrittr)
library(tidyverse)
library(reshape2)
library(S4Vectors)
library(SingleCellExperiment)
library(pheatmap)
library(png)
library(gridExtra)
library(knitr)
library(scales)
library(RColorBrewer)
library(Matrix.utils)
library(tibble)
library(ggplot2)
library(scater)
library(patchwork)
library(statmod)
library(ArchR)
library(clustree)
library(harmony)
library(gprofiler2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(AnnotationHub)
library(ReactomePA)
library(statmod)
library(edgeR)
library(speckle)
library(EnhancedVolcano)
library(decoupleR)
library(OmnipathR)
library(dorothea)
library(enrichplot)
library(png)
library(ezRun)
library(UpSetR)
library(ComplexHeatmap)

#Output paths
output_dir_data <- here::here("output", current_file)
if (!dir.exists(output_dir_data)) dir.create(output_dir_data)

if (!dir.exists(here::here("docs", "figure"))) dir.create(here::here("docs", "figure"))

output_dir_figs <- here::here("docs", "figure", paste0(current_file, ".Rmd"))
if (!dir.exists(output_dir_figs)) dir.create(output_dir_figs)

Load cardiomyocytes dataset

CM <- readRDS(here::here("output", "Cardiomyocytes_subclustering_26", "CM.rds"))
DefaultAssay(CM) <- "RNA"
CM <- NormalizeData(CM)

Differential abundance analysis

#Prepare matadata for differential abundance and differential expression testing
metadata <- CM@meta.data %>%
  dplyr::select(Sample:CO_l.min) %>%
  unique() %>%
  dplyr::select(-LVID_mm, -RVIT_mm) 
metadata <- metadata %>%
  dplyr::select(1:3, 5, 6, 4, 7:ncol(metadata)) %>%
  dplyr::arrange(match(Sample, names(table(CM$Sample))))

rownames(metadata) <- metadata$Sample

# Check if metadata$Sample matches names(table(CM$Sample))
all.equal(metadata$Sample, names(table(CM$Sample)))
[1] TRUE
ezInteractiveTableRmd(metadata)
# Filter for Group == "DCM"
metadata_DCM <- metadata %>%
  dplyr::filter(Group == "DCM")
# Initialize an empty data frame for storing differential abundance results
differential_abundance <- data.frame()
props <- getTransformedProps(CM$cell_state, CM$Sample, transform="logit")

# Loop through each metadata column from Age_years to CO_l.min
for (i in 6:(length(metadata_DCM))) {
  metadata_subset <- metadata_DCM[complete.cases(metadata_DCM[[i]]),]
  props_subset <- props
  props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
  props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
  props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]
  metadata_subset <- metadata_subset %>% 
    arrange(Sample, colnames(props_subset$Proportions))
  design <- model.matrix(~ metadata_subset[[i]])
  fit <- lmFit(props_subset$TransformedProps, design)
  fit <- eBayes(fit, robust=TRUE)
  differential_abundance_temp <- topTable(fit, n = Inf, coef = 2)
  differential_abundance_temp$metadata <- colnames(metadata_subset[i])
  differential_abundance_temp$cell_state <- rownames(differential_abundance_temp)
  differential_abundance <- rbind(differential_abundance, differential_abundance_temp)
}

# Clean up the environment
rm(metadata_subset, props_subset, design, fit, differential_abundance_temp)

# Reset row names and arrange the results by adjusted p-value
rownames(differential_abundance) <- NULL
differential_abundance <- dplyr::arrange(differential_abundance, adj.P.Val)

write.csv(differential_abundance, file = here::here(output_dir_data, "CM_Differential_Abundance.csv"), quote=F, row.names = F)

t <- ezInteractiveTableRmd(differential_abundance)
t[["x"]][["options"]][["pageLength"]] <- 10
t
metadata_subset <- metadata[complete.cases(metadata[["LVEDD_mm"]]),]

props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]

metadata_subset <- metadata_subset %>% 
  arrange(Sample, colnames(props_subset$Proportions))

design <- model.matrix(~ metadata_subset$LVEDD_mm)
fit.prop <- lmFit(props_subset$Proportions,design)

par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
  plot(metadata_subset$LVEDD_mm, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
       pch=16, cex=2, ylab="Proportions", xlab="LVEDD_mm", cex.lab=2, cex.axis=1.5,
       cex.main=2.5)
  abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
         lwd=1)
}

Version Author Date
61e8d17 GinoBonazza 2024-08-01
rm(metadata_subset, props_subset, design, fit.prop)
metadata_subset <- metadata[complete.cases(metadata[["LVEF_percent"]]),]

props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]

metadata_subset <- metadata_subset %>% 
  arrange(Sample, colnames(props_subset$Proportions))

design <- model.matrix(~ metadata_subset$LVEF_percent)
fit.prop <- lmFit(props_subset$Proportions,design)

par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
  plot(metadata_subset$LVEF_percent, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
       pch=16, cex=2, ylab="Proportions", xlab="LVEF_percent", cex.lab=2, cex.axis=1.5,
       cex.main=2.5)
  abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
         lwd=1)
}

Version Author Date
61e8d17 GinoBonazza 2024-08-01
rm(metadata_subset, props_subset, design, fit.prop)
metadata_subset <- metadata[complete.cases(metadata[["Age_years"]]),]

props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]

metadata_subset <- metadata_subset %>% 
  arrange(Sample, colnames(props_subset$Proportions))

design <- model.matrix(~ metadata_subset$Age_years)
fit.prop <- lmFit(props_subset$Proportions,design)

par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
  plot(metadata_subset$Age_years, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
       pch=16, cex=2, ylab="Proportions", xlab="Age_years", cex.lab=2, cex.axis=1.5,
       cex.main=2.5)
  abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
         lwd=1)
}

Version Author Date
61e8d17 GinoBonazza 2024-08-01
rm(metadata_subset, props_subset, design, fit.prop)
metadata_subset <- metadata[complete.cases(metadata[["RVSP_mmHg"]]),]

props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]

metadata_subset <- metadata_subset %>% 
  arrange(Sample, colnames(props_subset$Proportions))

design <- model.matrix(~ metadata_subset$RVSP_mmHg)
fit.prop <- lmFit(props_subset$Proportions,design)

par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
  plot(metadata_subset$RVSP_mmHg, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
       pch=16, cex=2, ylab="Proportions", xlab="RVSP_mmHg", cex.lab=2, cex.axis=1.5,
       cex.main=2.5)
  abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
         lwd=1)
}

Version Author Date
61e8d17 GinoBonazza 2024-08-01
rm(metadata_subset, props_subset, design, fit.prop)
design <- model.matrix(~ 0 + metadata$Group)
colnames(design) <- c("DCM", "HC")
mycontr <- makeContrasts(DCM-HC, levels=design)
differential_abundance_HC <- propeller.ttest(props, design, contrasts = mycontr, robust=TRUE, trend=FALSE, sort=TRUE)

differential_abundance_HC
    PropMean.DCM PropMean.HC PropRatio Tstatistic      P.Value         FDR
CM3   0.11135354  0.03836732 2.9023016   3.800582 0.0006739595 0.003369797
CM4   0.08662428  0.02831785 3.0589987   3.202724 0.0032613900 0.008153475
CM1   0.50911696  0.57790022 0.8809773  -1.527054 0.1374288591 0.229048098
CM2   0.24032230  0.29095689 0.8259722  -1.002228 0.3244063073 0.405507884
CM5   0.05258291  0.06445772 0.8157738  -0.479165 0.6353640356 0.635364036
par(mfrow = c(1, 5), mar = c(3, 5, 3, 1))
for (i in seq(1, 5, 1)) {
  stripchart(props$Proportions[i,] ~ metadata$Group,
             vertical = TRUE, pch = 16, method = "jitter", 
             col = c("orange", "purple"),
             cex = 2, ylab = "Proportions", cex.lab = 2, cex.axis = 2)
  title(rownames(props$Proportions)[i], cex.main = 2.5)
}

Version Author Date
61e8d17 GinoBonazza 2024-08-01
rm(design, mycontr)

Differential expression analysis

Remove genes expressed in <1% of cells and MT and RP genes

percent_stats <- Percent_Expressing(seurat_object = CM, features = rownames(CM), assay = "RNA", entire_object = TRUE)
percent_stats <- percent_stats %>%
  dplyr::rename(percentCells = All_Cells)
percent_stats$gene <- rownames(percent_stats)

write.csv(percent_stats, file = here::here(output_dir_data, "CM_percent_stats.csv"), quote=F, row.names = T)
percent_stats <- read.csv(here::here(output_dir_data, "CM_percent_stats.csv"))
rownames(percent_stats) <- percent_stats$gene

keep_genes <- rownames(dplyr::filter(percent_stats, percentCells > 1))
data <- CM[which(rownames(CM) %in% keep_genes),]
rm(keep_genes)

t <- ezInteractiveTableRmd(percent_stats)
t[["x"]][["options"]][["pageLength"]] <- 10
t

Create and prepare pseudobulk object for differential expression analysis

pseudocounts <- Seurat2PB(data, sample="Sample", cluster = "cell_type")
colnames(pseudocounts) <- pseudocounts$samples$sample
saveRDS(pseudocounts, 
        here::here(output_dir_data, "CM_pseudocounts.rds"))
keep_samples <- pseudocounts$samples$lib.size > 5e4
pseudocounts <- pseudocounts[, keep_samples]
keep_genes <- filterByExpr(pseudocounts)
pseudocounts <- pseudocounts[keep_genes, ]
rm(keep_samples, keep_genes, data)

Create empty objects for the DE output

results <- list()
signif <- list()
volcano <- list()
n_de_genes <- data.frame()
pseudocounts_DCM <- pseudocounts[, metadata_DCM$Sample]
all.equal(metadata_DCM$Sample, colnames(pseudocounts_DCM$counts))
[1] TRUE
for (i in 6:(length(metadata_DCM))) {
  metadata_subset <- metadata_DCM[complete.cases(metadata_DCM[[i]]),]
  metadata_subset[,i] <- rescale(metadata_subset[,i])
  design <- model.matrix(~ metadata_subset[[i]])
  colnames(design)[2] <- colnames(metadata_subset[i])
  
  pseudocounts_subset <- pseudocounts_DCM[, metadata_subset$Sample]
  pseudocounts_subset <- normLibSizes(pseudocounts_subset)
  pseudocounts_subset <- estimateDisp(pseudocounts_subset, design)
  
  fit <- glmQLFit(pseudocounts_subset, design, robust=TRUE)
  fit <- glmQLFTest(fit, coef = 2)
  
  results[[i-5]] <- topTags(fit, n = Inf)$table
  results[[i-5]] <- merge(results[[i-5]], percent_stats[, c("gene", "percentCells")], by = "gene", all.x = FALSE)
  signif[[i-5]] <- results[[i-5]] %>% dplyr::filter(FDR < 0.05, abs(logFC) > 0.5) %>% 
    dplyr::arrange(FDR)

  names(results)[i-5] <- colnames(metadata_subset[i])
  names(signif)[i-5] <- colnames(metadata_subset[i])
  
  write.csv(results[[i-5]], file = here::here(output_dir_data, paste0("CM_DE_Results_", colnames(metadata[i]), ".csv")), quote=F, row.names = F)
  write.csv(signif[[i-5]], file = here::here(output_dir_data, paste0("CM_DE_Significant_", colnames(metadata[i]), ".csv")), quote=F, row.names = F)
  
  n_de_genes_temp <- data.frame(metadata = colnames(metadata_subset[i]),
                                n_upregulated = sum(signif[[i-5]]$logFC > 0.5),
                                n_downregulated = sum(signif[[i-5]]$logFC < -0.5),
                                n_tot = nrow(signif[[i-5]])
                                )
  n_de_genes <- rbind(n_de_genes, n_de_genes_temp)
  
  volcano[[i-5]] <- EnhancedVolcano(results[[i-5]],
                  lab = results[[i-5]]$gene,
                  x = "logFC",
                  y = "FDR",
                  labSize = 0,
                  titleLabSize = 16,
                  subtitleLabSize = 14,
                  axisLabSize = 12,
                  captionLabSize = 9,
                  pointSize = 0.5,
                  FCcutoff = 0.5,
                  pCutoff  = 0.05,
                  ylim = c(0, 3.5),
                  col = c("black", "black", "black", "red"),
                  colAlpha = 1,
                  drawConnectors = FALSE,
                  subtitle = NULL,
                  title = paste0("Cardiomyocytes", "\n", colnames(metadata_subset[i]), "\nn = ", nrow(metadata_subset))
  ) + theme(legend.position = "none")
  names(volcano)[i-5] <- colnames(metadata_subset[i])
}
rm(metadata_subset, design, pseudocounts_subset, fit, n_de_genes_temp)
print((volcano[[1]] | volcano[[2]] | volcano[[3]] | volcano[[4]] | volcano[[5]] | volcano[[6]] | volcano[[7]]) / 
  (volcano[[8]] | volcano[[9]] | volcano[[10]] |  volcano[[11]] | volcano[[12]] | volcano[[13]] | volcano[[14]]))

Version Author Date
61e8d17 GinoBonazza 2024-08-01
gene_sets <- list()
for (i in 1:length(signif)) {
  gene_sets[[i]] <- signif[[i]]$gene
  names(gene_sets)[i] <- names(signif)[i]
  }
m <- make_comb_mat(gene_sets, mode = "intersect", min_set_size = 10)
m <- m[comb_degree(m) >= 2]
m <- m[comb_size(m) >= 10]
upset_plot <- UpSet(m, column_title = "Cardiomyocytes", comb_order = order(comb_size(m)), 
                    top_annotation = upset_top_annotation(m, add_numbers = TRUE),
                    right_annotation = upset_right_annotation(m, add_numbers = TRUE))
draw(upset_plot)

Over-representation analysis

# Calculate gene counts across cells
gene_counts <- rowSums(CM@assays$RNA@counts > 0)

# Filter genes to include only those expressed in at least 3 cells
universe_genes <- names(gene_counts[gene_counts >= 3])
write.csv(as.data.frame(universe_genes), here::here(output_dir_data, "CM_universe_genes.csv"), quote=F, row.names = F)

# Convert universe gene symbols to Entrez IDs
universe_entrez <- bitr(universe_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID
write.csv(as.data.frame(universe_entrez), here::here(output_dir_data, "CM_universe_entrez.csv"), quote=F, row.names = F)

rm(gene_counts)
up_genes <- filter(signif[["mPAP_mmHg"]], logFC > 0.5)$gene
up_genes <- bitr(up_genes, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
down_genes <- filter(signif[["mPAP_mmHg"]], logFC < -0.5)$gene
down_genes <- bitr(down_genes, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
reorder_GO_by_pvalue <- function(enrichGO_result) {
  go_results <- enrichGO_result@result
  go_results_sorted <- go_results[order(go_results$p.adjust), ]
  enrichGO_result@result <- go_results_sorted
  return(enrichGO_result)
}
mPAP_GO_up <- enrichGO(up_genes, OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
mPAP_GO_up <- gsfilter(mPAP_GO_up, by = 'Count', min = 3)
mPAP_GO_up <- reorder_GO_by_pvalue(mPAP_GO_up)
saveRDS(mPAP_GO_up, here::here(output_dir_data, "CM_ORA_GO_mPAP_up.rds"))

mPAP_GO_down <- enrichGO(down_genes, OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
mPAP_GO_down <- gsfilter(mPAP_GO_down, by = 'Count', min = 3)
mPAP_GO_down <- reorder_GO_by_pvalue(mPAP_GO_down)
saveRDS(mPAP_GO_down, here::here(output_dir_data, "CM_ORA_GO_mPAP_down.rds"))
mPAP_GO_up <- readRDS(here::here(output_dir_data, "CM_ORA_GO_mPAP_up.rds"))
mPAP_GO_down <- readRDS(here::here(output_dir_data, "CM_ORA_GO_mPAP_down.rds"))
p1 <- dotplot(mPAP_GO_up, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
p2 <- dotplot(mPAP_GO_down, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nLow mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
print(p1|p2)

mPAP_GO_simplified_up <- simplify(
  mPAP_GO_up,
  cutoff = 0.75,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
mPAP_GO_simplified_down <- simplify(
  mPAP_GO_down,
  cutoff = 0.75,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
p1 <- dotplot(mPAP_GO_simplified_up, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
p2 <- dotplot(mPAP_GO_simplified_down, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nLow mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
print(p1 | p2)

mPAP_KEGG_up <- enrichKEGG(up_genes, organism = 'hsa', universe = universe_entrez)
mPAP_KEGG_up <- gsfilter(mPAP_KEGG_up, by = 'Count', min = 3)
saveRDS(mPAP_KEGG_up, here::here(output_dir_data, "CM_ORA_KEGG_mPAP_up.rds"))
mPAP_KEGG_down <- enrichKEGG(down_genes,  organism = 'hsa', universe = universe_entrez)
mPAP_KEGG_down <- gsfilter(mPAP_KEGG_down, by = 'Count', min = 3)
mPAP_KEGG_up <- readRDS(here::here(output_dir_data, "CM_ORA_KEGG_mPAP_up.rds"))
p1 <- dotplot(mPAP_KEGG_up, showCategory = 10, title = paste0("KEGG Pathways", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
print(p1)

mPAP_REACTOME_up <- enrichPathway(up_genes, organism = 'human', universe = universe_entrez, readable = TRUE, pvalueCutoff = 0.1)
mPAP_REACTOME_up <- gsfilter(mPAP_REACTOME_up, by = 'Count', min = 3)
saveRDS(mPAP_REACTOME_up, here::here(output_dir_data, "CM_ORA_REACTOME_mPAP_up.rds"))
mPAP_REACTOME_down <- enrichPathway(down_genes, organism = 'human', universe = universe_entrez, readable = TRUE, pvalueCutoff = 0.1)
mPAP_REACTOME_down <- gsfilter(mPAP_REACTOME_down, by = 'Count', min = 3)
mPAP_REACTOME_up <- readRDS(here::here(output_dir_data, "CM_ORA_REACTOME_mPAP_up.rds"))
p1 <- dotplot(mPAP_REACTOME_up, showCategory = 10, title = paste0("REACTOME", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))
print(p1)

rm(mPAP_GO_up, mPAP_GO_down, mPAP_GO_simplified_up, mPAP_GO_simplified_down, mPAP_KEGG_up, mPAP_KEGG_down, mPAP_REACTOME_up, mPAP_REACTOME_down)

Gene set enrichment analysis

ranked_genes <- results[["mPAP_mmHg"]]
ranked_genes$metric <- ranked_genes$logFC*-log10(ranked_genes$PValue)
entrezid <- bitr(ranked_genes$gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db", drop = TRUE)
ranked_genes <- merge(ranked_genes, entrezid, by.x = "gene", by.y = "SYMBOL")
ranked_genes <- ranked_genes[!duplicated(ranked_genes$metric), ]
genelist <- ranked_genes$metric
names(genelist) <- ranked_genes$ENTREZID
genelist <- genelist[order(genelist, decreasing = T)]
write.csv(ranked_genes, file = here::here(output_dir_data, "CM_ranked_genes.csv"), quote=F, row.names = F)
t <- ezInteractiveTableRmd(ranked_genes)
t[["x"]][["options"]][["pageLength"]] <- 10
t
gseaGO_all <- gseGO(geneList = genelist,
                OrgDb = org.Hs.eg.db,
                ont = "ALL",
                keyType = "ENTREZID",
                nPermSimple = 1000000,
                minGSSize = 10)
gseaGO_all <- reorder_GO_by_pvalue(gseaGO_all)

saveRDS(gseaGO_all, here::here(output_dir_data, "CM_GSEA_GO_mPAP.rds"))
gseaGO_all <- readRDS(here::here(output_dir_data, "CM_GSEA_GO_mPAP.rds"))
gseaGO_all_simplified <- simplify(
  gseaGO_all,
  cutoff = 0.75,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
p1 <- dotplot(gseaGO_all_simplified, showCategory = 10, split=".sign", title = paste0("GO Enrichment Analysis", "\n mPAP-associated genes"), label_format = 30, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"),
        strip.text = element_text(face = "bold", size = 18)) +
  facet_grid(.~.sign)
print(p1)

classify_terms <- function(df) {
  df$sign <- ifelse(df$NES > 0, "activated", "suppressed")
  return(df)
}
gseaGO_all_simplified@result <- classify_terms(gseaGO_all_simplified@result)
# Split the GSEA results
activated_terms <- subset(gseaGO_all_simplified@result, sign == "activated")
suppressed_terms <- subset(gseaGO_all_simplified@result, sign == "suppressed")

# Create two separate enrichResult objects
gseaGO_all_activated <- gseaGO_all_simplified
gseaGO_all_activated@result <- activated_terms
gseaGO_all_activated <- pairwise_termsim(gseaGO_all_activated)

gseaGO_all_suppressed <- gseaGO_all_simplified
gseaGO_all_suppressed@result <- suppressed_terms
gseaGO_all_suppressed <- pairwise_termsim(gseaGO_all_suppressed)

# Plot activated terms
p1 <- emapplot(gseaGO_all_activated, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("Activated Terms") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Plot suppressed terms
p2 <- emapplot(gseaGO_all_suppressed, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("Suppressed Terms") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Combine the plots
combined_plot <- p1 | p2
print(combined_plot)

p1 <- ridgeplot(gseaGO_all_simplified, showCategory = 20) + labs(x = "enrichment distribution")
print(p1)

gseaKEGG <- gseKEGG(geneList = genelist,
                organism = "hsa",
                keyType = "kegg",
                nPermSimple = 100000)
saveRDS(gseaKEGG, here::here(output_dir_data, "CM_GSEA_KEGG_mPAP.rds"))
gseaKEGG <- readRDS(here::here(output_dir_data, "CM_GSEA_KEGG_mPAP.rds"))
p1 <- dotplot(gseaKEGG, showCategory = 10, split=".sign", title = paste0("GSEA\nKEGG Pathways", "\n mPAP-associated genes"), label_format = 30, font.size = 15) +
  theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"),
        strip.text = element_text(face = "bold", size = 18)) +
  facet_grid(.~.sign)
print(p1)

gseaKEGG@result <- classify_terms(gseaKEGG@result)
# Split the GSEA results
activated_terms <- subset(gseaKEGG@result, sign == "activated")
suppressed_terms <- subset(gseaKEGG@result, sign == "suppressed")

# Create two separate enrichResult objects
gseaKEGG_activated <- gseaKEGG
gseaKEGG_activated@result <- activated_terms

gseaKEGG_suppressed <- gseaKEGG
gseaKEGG_suppressed@result <- suppressed_terms

# Plot activated terms
p1 <- emapplot(pairwise_termsim(gseaKEGG_activated), showCategory = 20, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("Activated Terms") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Plot suppressed terms
p2 <- emapplot(pairwise_termsim(gseaKEGG_suppressed), showCategory = 20, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("Suppressed Terms") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Combine the plots
combined_plot <- p1 | p2
print(combined_plot)

p1 <- ridgeplot(gseaKEGG, showCategory = 20) + labs(x = "enrichment distribution")
print(p1)

gseaGO_all_activated_df <- setReadable(gseaGO_all_activated, OrgDb = org.Hs.eg.db)
gseaGO_all_activated_df <- gseaGO_all_activated_df@result

Heatmaps

metadata_mPAP <- metadata[complete.cases(metadata[["mPAP_mmHg"]]),]

pseudocounts_mPAP <- pseudocounts[, metadata_mPAP$Sample]
pseudocounts_mPAP <- normLibSizes(pseudocounts_mPAP)

counts_mPAP <- cpm(pseudocounts_mPAP, log = TRUE)
## Set a color-blind friendly palette
heat_colors <- rev(brewer.pal(11, "PuOr"))

metadata_heatmap <- metadata_mPAP %>% dplyr::select(Sample, mPAP_mmHg, Batch_processing, Age_years) %>% dplyr::arrange(mPAP_mmHg)

DEGs_mPAP <- signif[["mPAP_mmHg"]]

sig_counts <- counts_mPAP[rownames(counts_mPAP) %in% DEGs_mPAP$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = TRUE, 
         show_rownames = FALSE,
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 10, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 10, 
         height = 20) 

Version Author Date
61e8d17 GinoBonazza 2024-08-01
rm(metadata_heatmap, sig_counts)
metadata_heatmap <- metadata_mPAP %>% dplyr::select(Sample, mPAP_mmHg, Batch_processing, Age_years) %>% dplyr::arrange(mPAP_mmHg)

DEGs_mPAP_top <- DEGs_mPAP
DEGs_mPAP_top <- DEGs_mPAP_top %>% 
  dplyr::filter(percentCells > 5) 
DEGs_mPAP_top_up <- DEGs_mPAP_top %>% 
  dplyr::filter(logFC > 2) %>% 
  dplyr::arrange(FDR)
DEGs_mPAP_top_down <- DEGs_mPAP_top %>% 
  dplyr::filter(logFC < -2) %>% 
  dplyr::arrange(FDR)

DEGs_mPAP_top <- dplyr::bind_rows(head(DEGs_mPAP_top_up, 20), head(DEGs_mPAP_top_down, 20))
rm(DEGs_mPAP_top_up, DEGs_mPAP_top_down)

sig_counts <- counts_mPAP[rownames(counts_mPAP) %in% DEGs_mPAP_top$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = TRUE, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 10, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 15, 
         height = 20,
         show_rownames = T) 

Version Author Date
61e8d17 GinoBonazza 2024-08-01
rm(metadata_heatmap, sig_counts)
metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Batch_processing, Age_years)

pseudocounts_mPAP_HC <- pseudocounts[, metadata_heatmap$Sample]
pseudocounts_mPAP_HC <- normLibSizes(pseudocounts_mPAP_HC)

counts_mPAP_HC <- cpm(pseudocounts_mPAP_HC, log = TRUE)

sig_counts <- counts_mPAP_HC[rownames(counts_mPAP_HC) %in% DEGs_mPAP$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = TRUE, 
         show_rownames = FALSE,
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 10, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 10, 
         height = 20) 

Version Author Date
61e8d17 GinoBonazza 2024-08-01
rm(metadata_heatmap, sig_counts)
metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Batch_processing, Age_years)

pseudocounts_mPAP_HC <- pseudocounts[, metadata_heatmap$Sample]
pseudocounts_mPAP_HC <- normLibSizes(pseudocounts_mPAP_HC)

counts_mPAP_HC <- cpm(pseudocounts_mPAP_HC, log = TRUE)

sig_counts <- counts_mPAP_HC[rownames(counts_mPAP_HC) %in% DEGs_mPAP_top$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = TRUE, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 10, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 15, 
         height = 20,
         show_rownames = T) 

Version Author Date
61e8d17 GinoBonazza 2024-08-01

DEGs hierarchical clustering

# Filter and arrange metadata
metadata_mPAP_HC <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC") %>%
  dplyr::mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  dplyr::arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Batch_processing, Age_years)

# Normalize pseudocounts and calculate log CPM
pseudocounts_mPAP_HC <- pseudocounts[, metadata_mPAP_HC$Sample]
pseudocounts_mPAP_HC <- normLibSizes(pseudocounts_mPAP_HC)
counts_mPAP_HC <- edgeR::cpm(pseudocounts_mPAP_HC, log = TRUE)

# Filter counts matrix by DEGs
counts_matrix <- counts_mPAP_HC[rownames(counts_mPAP_HC) %in% DEGs_mPAP$gene, ]

# Convert to long format and add Condition dividing DCM samples in 3 groups based on mPAP
counts_tibble <- counts_matrix %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "gene") %>%
  tidyr::pivot_longer(-gene, names_to = "sample", values_to = "log_pseudocounts") %>%
  dplyr::mutate(Condition = dplyr::case_when(
    sample %in% c("CZD4", "CZD5", "CZD6", "CZD7") ~ "HC",
    sample %in% c("PK439", "PK375", "PK229", "P5", "P8", "PK153") ~ "DCM - Mild PH",
    sample %in% c("PK385", "PK437", "PK323", "PK403", "PK213", "PK461") ~ "DCM - Moderate PH",
    sample %in% c("PK369", "PK357", "PK337", "PK345", "PK401", "PK463") ~ "DCM - Severe PH"
  ))

# Compute the mean log counts for each gene and condition
mean_tibble <- counts_tibble %>%
  dplyr::group_by(gene, Condition) %>%
  dplyr::summarise(mean_log_counts = mean(log_pseudocounts, na.rm = TRUE), .groups = 'drop')

# Convert to wide format
wide_mean_tibble <- mean_tibble %>%
  tidyr::pivot_wider(names_from = Condition, values_from = mean_log_counts)

# Convert the wide tibble to a matrix and transpose
mean_matrix <- wide_mean_tibble %>%
  tibble::column_to_rownames(var = "gene") %>%
  as.matrix()

mean_matrix  <- mean_matrix %>%
  t() %>% 
  scale() %>% 
  t()

# Perform hierarchical clustering
gene_dist <- dist(mean_matrix)
gene_hclust <- hclust(gene_dist, method = "complete")

# Plot the dendrogram
plot(gene_hclust, labels = FALSE)
abline(h = 3, col = "brown", lwd = 2) # Add horizontal line to illustrate cutting dendrogram

Version Author Date
61e8d17 GinoBonazza 2024-08-01
# Cut the dendrogram into 4 clusters
gene_cluster <- cutree(gene_hclust, k = 4) %>%
  tibble::enframe() %>%
  dplyr::rename(gene = name, cluster = value)

# Convert scaled matrix to long format and join with gene_cluster
mean_tibble_cluster <- mean_matrix %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var = "gene") %>%
  tidyr::pivot_longer(-gene, names_to = "Condition", values_to = "scaled_log_pseudocounts") %>%
  dplyr::inner_join(gene_cluster, by = "gene")

# Set factor levels for Condition
mean_tibble_cluster$Condition <- factor(mean_tibble_cluster$Condition, levels = c("HC", "DCM - Mild PH", "DCM - Moderate PH", "DCM - Severe PH"))

# Calculate the number of genes per cluster
cluster_counts <- mean_tibble_cluster %>%
  dplyr::group_by(cluster) %>%
  dplyr::summarise(gene_count = dplyr::n_distinct(gene))

# Create a named vector for custom facet labels
custom_labels <- setNames(paste("Cluster", cluster_counts$cluster, "-", cluster_counts$gene_count, "genes"), cluster_counts$cluster)

# Plot the data with custom facet labels
mean_tibble_cluster %>%
  ggplot(aes(Condition, scaled_log_pseudocounts)) +
  geom_line(aes(group = gene), alpha = 0.2) +
  geom_line(stat = "summary", fun = "median", colour = "brown", size = 1.5, aes(group = 1)) +
  facet_grid(cols = vars(cluster), labeller = as_labeller(custom_labels)) +
  theme(
    strip.text = element_text(face = "bold", size = 12),       # Make facet labels bold
    axis.text.x = element_text(angle = 45, hjust = 1, size = 10) # Rotate x-axis labels by 45 degrees
  )

Version Author Date
61e8d17 GinoBonazza 2024-08-01

Over-representation DEGs clusters

DEGs_clusters <- list()
DEGs_clusters[["cluster1"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "1")$gene)
DEGs_clusters[["cluster2"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "2")$gene)
DEGs_clusters[["cluster3"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "3")$gene)
DEGs_clusters[["cluster4"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "4")$gene)
DEGs_clusters_entrez <- list()
DEGs_clusters_entrez[["cluster1"]] <- bitr(DEGs_clusters[["cluster1"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
DEGs_clusters_entrez[["cluster2"]] <- bitr(DEGs_clusters[["cluster2"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
DEGs_clusters_entrez[["cluster3"]] <- bitr(DEGs_clusters[["cluster3"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
DEGs_clusters_entrez[["cluster4"]] <- bitr(DEGs_clusters[["cluster4"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
cluster1_GO <- enrichGO(DEGs_clusters_entrez[["cluster1"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster1_GO <- gsfilter(cluster1_GO, by = 'Count', min = 3)
nrow(cluster1_GO)
cluster2_GO <- enrichGO(DEGs_clusters_entrez[["cluster2"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster2_GO <- gsfilter(cluster2_GO, by = 'Count', min = 3)
nrow(cluster2_GO)
cluster3_GO <- enrichGO(DEGs_clusters_entrez[["cluster3"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster3_GO <- gsfilter(cluster3_GO, by = 'Count', min = 10)
nrow(cluster3_GO)
cluster3_GO <- simplify(
  cluster3_GO,
  cutoff = 0.7,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
cluster3_GO <- reorder_GO_by_pvalue(cluster3_GO)
saveRDS(cluster3_GO, here::here(output_dir_data, "CM_cluster3_ORA.rds"))
cluster3_GO <- readRDS(here::here(output_dir_data, "CM_cluster3_ORA.rds"))
dotplot(cluster3_GO, 
        showCategory = 15, 
        title = paste0("GO - Over-representation analysis", "\nDEGs cluster 3"), 
        label_format = 27, 
        font.size = 15) +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))

cluster4_GO <- enrichGO(DEGs_clusters_entrez[["cluster4"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster4_GO <- gsfilter(cluster4_GO, by = 'Count', min = 10)
nrow(cluster4_GO)
cluster4_GO <- simplify(
  cluster4_GO,
  cutoff = 0.7,
  by = "p.adjust",
  select_fun = min,
  measure = "Wang",
  semData = NULL
)
cluster4_GO <- reorder_GO_by_pvalue(cluster4_GO)
saveRDS(cluster4_GO, here::here(output_dir_data, "CM_cluster4_ORA.rds"))
cluster4_GO <- readRDS(here::here(output_dir_data, "CM_cluster4_ORA.rds"))
dotplot(cluster4_GO, 
        showCategory = 15, 
        title = paste0("GO - Over-representation analysis", "\nDEGs cluster 4"), 
        label_format = 27, 
        font.size = 15) +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5),
        axis.text.y = element_text(face = "bold"))

cluster3_GO <- pairwise_termsim(cluster3_GO)
cluster4_GO <- pairwise_termsim(cluster4_GO)

p1 <- emapplot(cluster3_GO, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("DEGs cluster 3") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Plot suppressed terms
p2 <- emapplot(cluster4_GO, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
  ggtitle("DEGs cluster 4") +
  theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))

# Combine the plots
combined_plot <- p1 | p2
print(combined_plot)

Top genes

DEGs_mPAP_cluster1 <- DEGs_mPAP %>% 
  filter(gene %in% DEGs_clusters[["cluster1"]])
top_DEGs_mPAP_cluster1 <- DEGs_mPAP_cluster1 %>% 
  filter(percentCells > 5) %>%
  filter(FDR < 0.01) %>%
  arrange(logFC) %>%
  head(20)

DEGs_mPAP_cluster2 <- DEGs_mPAP %>% 
  filter(gene %in% DEGs_clusters[["cluster2"]])
top_DEGs_mPAP_cluster2 <- DEGs_mPAP_cluster2 %>% 
  filter(percentCells > 5) %>%
  filter(FDR < 0.01) %>%
  arrange(logFC) %>%
  head(20)

DEGs_mPAP_cluster3 <- DEGs_mPAP %>% 
  filter(gene %in% DEGs_clusters[["cluster3"]])
top_DEGs_mPAP_cluster3 <- DEGs_mPAP_cluster3 %>% 
  filter(percentCells > 5) %>%
  filter(FDR < 0.01) %>%
  arrange(desc(logFC)) %>%
  head(20)

DEGs_mPAP_cluster4 <- DEGs_mPAP %>% 
  filter(gene %in% DEGs_clusters[["cluster4"]])
top_DEGs_mPAP_cluster4 <- DEGs_mPAP_cluster4 %>% 
  filter(percentCells > 5) %>%
  filter(FDR < 0.01) %>%
  arrange(desc(logFC)) %>%
  head(20)
metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Age_years)

sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster1$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = F, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 7, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 11, 
         cellheight = 17,
         show_rownames = T,
         main = "Top cluster 1 DEGs") 

Version Author Date
61e8d17 GinoBonazza 2024-08-01
ezInteractiveTableRmd(top_DEGs_mPAP_cluster1)
metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Age_years)

sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster2$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = F, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 7, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 11, 
         cellheight = 17,
         show_rownames = T,
         main = "Top cluster 2 DEGs") 

Version Author Date
61e8d17 GinoBonazza 2024-08-01
ezInteractiveTableRmd(top_DEGs_mPAP_cluster2)
metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Age_years)

sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster3$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = F, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 7, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 11, 
         cellheight = 17,
         show_rownames = T,
         main = "Top cluster 3 DEGs") 

Version Author Date
61e8d17 GinoBonazza 2024-08-01
ezInteractiveTableRmd(top_DEGs_mPAP_cluster3)
metadata_heatmap <- metadata %>% 
  dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>% 
  mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
  arrange(Group, mPAP_mmHg) %>%
  dplyr::select(Sample, Group, mPAP_mmHg, Age_years)

sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster4$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]

## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts, 
         color = heat_colors, 
         cluster_rows = F, 
         annotation_col = metadata_heatmap, 
         border_color = NA, 
         fontsize = 7, 
         scale = "row",
         cluster_cols=F,
         fontsize_row = 11, 
         cellheight = 17,
         show_rownames = T,
         main = "Top cluster 4 DEGs") 

Version Author Date
61e8d17 GinoBonazza 2024-08-01
ezInteractiveTableRmd(top_DEGs_mPAP_cluster4)
signif_cluster <- signif[["mPAP_mmHg"]] %>%
  left_join(gene_cluster) %>%
  left_join(ranked_genes)
write.csv(signif_cluster, file = here::here(output_dir_data, "CM_signif_cluster.csv"), quote=F, row.names = F)
ezInteractiveTableRmd(signif_cluster)

Scores

terms <- list(
  muscle_structure_development = c("actin", "muscle", "sarcomere", "cardiac", "fibril", "fiber", "actomyosin"),
  autophagy = c("autophag", "disassembly", "catabolic", "starvation"),
  transcription = c("transcription", "receptor"),
  vesicular_transport = c("transport", "Golgi", "vacuo", "vesicle", "endosom")
)

get_unique_genes <- function(df, terms) {
  filtered_df <- df %>%
    filter(grepl(paste(terms, collapse = "|"), Description, ignore.case = TRUE))
  unique(unlist(strsplit(filtered_df$core_enrichment, split = "/")))
}

genes_of_interest <- lapply(terms, get_unique_genes, df = gseaGO_all_activated_df)
names(genes_of_interest) <- names(terms)
saveRDS(genes_of_interest, here::here(output_dir_data, "CM_genes_of_interest.rds"))
for (i in 1:length(genes_of_interest)) {
  score_genes <- list(genes_of_interest[[i]])
  score_name <- paste0(names(genes_of_interest)[i], "_score")
  CM <- AddModuleScore(object = CM, features = score_genes, name = score_name)
}
metadata_CM <- CM@meta.data
colnames(metadata_CM) <- gsub("_score1", "_score", colnames(metadata_CM))
CM@meta.data <- metadata_CM
rm(metadata_CM)
plot_list <- lapply(paste0(names(genes_of_interest), "_score"), function(feature) {
  p <- FeaturePlot(CM, features = feature, order = FALSE)
  p + scale_color_viridis(discrete = FALSE, option = "turbo")
})

CombinePlots(plots = plot_list, ncol = 4)

Version Author Date
61e8d17 GinoBonazza 2024-08-01
plot_list <- lapply(paste0(names(genes_of_interest), "_score"), function(feature) {
  p <- FeaturePlot(CM, features = feature, order = TRUE)
  p + scale_color_viridis(discrete = FALSE, option = "turbo")
})

CombinePlots(plots = plot_list, ncol = 4)

Version Author Date
61e8d17 GinoBonazza 2024-08-01
VlnPlot(CM, features = paste0(names(genes_of_interest), "_score"), ncol = 4, group.by = "cell_state", pt.size = 0)

Version Author Date
61e8d17 GinoBonazza 2024-08-01
CM@meta.data[CM@meta.data$Sample %in% c("CZD4", "CZD5", "CZD6", "CZD7"), "Condition"] = "HC"
CM@meta.data[CM@meta.data$Sample %in% c("PK439", "PK375", "PK229", "P5", "P8", "PK153"), "Condition"] = "DCM - Mild PH"
CM@meta.data[CM@meta.data$Sample %in% c("PK385", "PK437", "PK323", "PK403", "PK213", "PK461"), "Condition"] = "DCM - Moderate PH"
CM@meta.data[CM@meta.data$Sample %in% c("PK369", "PK357", "PK337", "PK345", "PK401", "PK463"), "Condition"] = "DCM - Severe PH"
CM@meta.data[CM@meta.data$Sample %in% c("P3", "PK415", "PK441"), "Condition"] = "DCM - no mPAP data"

CM$Condition <- factor(CM$Condition, levels = c("HC", "DCM - Mild PH", "DCM - Moderate PH", "DCM - Severe PH", "DCM - no mPAP data"))
VlnPlot(CM, features = paste0(names(genes_of_interest), "_score"), ncol = 4, group.by = "Condition", pt.size = 0)

Version Author Date
61e8d17 GinoBonazza 2024-08-01

Transcription factor inference analysis

net <- get_collectri(organism='human', split_complexes=FALSE)
net
# A tibble: 43,178 × 3
   source target   mor
   <chr>  <chr>  <dbl>
 1 MYC    TERT       1
 2 SPI1   BGLAP      1
 3 SMAD3  JUN        1
 4 SMAD4  JUN        1
 5 STAT5A IL2        1
 6 STAT5B IL2        1
 7 RELA   FAS        1
 8 WT1    NR0B1      1
 9 NR0B2  CASP1      1
10 SP1    ALDOA      1
# ℹ 43,168 more rows
ranked_genes_TF <- ranked_genes %>%
  dplyr::select(gene, logFC, metric, PValue) %>% 
  remove_rownames() %>%
  column_to_rownames(var = "gene") %>%
  as.matrix()
head(ranked_genes_TF)
               logFC      metric       PValue
A1CF      -4.4846666 -20.7823457 2.322259e-05
A2M       -1.6943794  -6.0929917 2.535118e-04
A2M-AS1    0.2481086   0.1532492 2.411737e-01
A2ML1     -0.4535097  -0.1800881 4.007756e-01
A2ML1-AS1 -2.5607983  -8.7931153 3.683495e-04
A4GALT     2.0319539   6.1047409 9.899886e-04
# Run ulm
contrast_acts <- run_ulm(mat=ranked_genes_TF[, 'metric', drop=FALSE], net=net, .source='source', .target='target',
                  .mor='mor', minsize = 5)
head(contrast_acts)
# A tibble: 6 × 5
  statistic source condition  score p_value
  <chr>     <chr>  <chr>      <dbl>   <dbl>
1 ulm       ABL1   metric     1.49    0.137
2 ulm       AHR    metric     1.06    0.287
3 ulm       AIP    metric     0.893   0.372
4 ulm       AIRE   metric    -0.730   0.465
5 ulm       AP1    metric     0.776   0.438
6 ulm       APEX1  metric     0.458   0.647
n_tfs <- 40

# Filter top TFs in both signs
f_contrast_acts <- contrast_acts %>%
  mutate(rnk = NA)
msk <- f_contrast_acts$score > 0
f_contrast_acts[msk, 'rnk'] <- rank(-f_contrast_acts[msk, 'score'])
f_contrast_acts[!msk, 'rnk'] <- rank(-abs(f_contrast_acts[!msk, 'score']))
tfs <- f_contrast_acts %>%
  arrange(rnk) %>%
  head(n_tfs) %>%
  pull(source)
f_contrast_acts <- f_contrast_acts %>%
  filter(source %in% tfs)

# Plot
ggplot(f_contrast_acts, aes(x = reorder(source, score), y = score)) + 
    geom_bar(aes(fill = score), stat = "identity") +
    scale_fill_gradient2(low = "darkblue", high = "indianred", 
        mid = "whitesmoke", midpoint = 0) + 
    theme_minimal() +
    theme(axis.title = element_text(face = "bold", size = 12),
        axis.text.x = 
            element_text(angle = 45, hjust = 1, size =10, face= "bold"),
        axis.text.y = element_text(size =10, face= "bold"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
    xlab("TFs")

Version Author Date
61e8d17 GinoBonazza 2024-08-01
FeaturePlot(CM, features = f_contrast_acts$source, ncol = 8)

Version Author Date
61e8d17 GinoBonazza 2024-08-01
contrast_acts_filtered <- contrast_acts %>%
  filter(source %in% filter(percent_stats, percentCells > 5)$gene)
n_tfs <- 40

# Filter top TFs in both signs
f_contrast_acts <- contrast_acts_filtered %>%
  mutate(rnk = NA) %>%
  filter(p_value < 0.05)
msk <- f_contrast_acts$score > 0
f_contrast_acts[msk, 'rnk'] <- rank(-f_contrast_acts[msk, 'score'])
f_contrast_acts[!msk, 'rnk'] <- rank(-abs(f_contrast_acts[!msk, 'score']))
tfs <- f_contrast_acts %>%
  arrange(rnk) %>%
  head(n_tfs) %>%
  pull(source)
f_contrast_acts <- f_contrast_acts %>%
  filter(source %in% tfs) 

# Plot
ggplot(f_contrast_acts, aes(x = reorder(source, score), y = score)) + 
    geom_bar(aes(fill = score), stat = "identity") +
    scale_fill_gradient2(low = "darkblue", high = "indianred", 
        mid = "whitesmoke", midpoint = 0) + 
    theme_minimal() +
    theme(axis.title = element_text(face = "bold", size = 12),
        axis.text.x = 
            element_text(angle = 45, hjust = 1, size =10, face= "bold"),
        axis.text.y = element_text(size =10, face= "bold"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
    xlab("TFs")

Version Author Date
61e8d17 GinoBonazza 2024-08-01
FeaturePlot(CM, features = f_contrast_acts$source, ncol = 6)

Version Author Date
61e8d17 GinoBonazza 2024-08-01
tf <- 'MYOCD'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

Version Author Date
61e8d17 GinoBonazza 2024-08-01
tf <- 'SRF'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

tf <- 'MEF2A'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

Version Author Date
182f9de GinoBonazza 2024-08-19
5cd0f74 GinoBonazza 2024-08-01
61e8d17 GinoBonazza 2024-08-01

https://www.nature.com/articles/s41419-023-05665-8

tf <- 'MEF2C'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

Version Author Date
182f9de GinoBonazza 2024-08-19
61e8d17 GinoBonazza 2024-08-01
tf <- 'NFATC3'

df <- net %>%
  filter(source == tf) %>%
  arrange(target) %>%
  mutate(ID = target, color = "3") %>%
  column_to_rownames('target')

inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
  mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
  mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
  mutate(color = if_else(mor < 0 & t_value < 0, '1', color))

ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
  geom_point() +
  scale_colour_manual(values = c("red","royalblue3","grey")) +
  geom_label_repel(aes(label = ID, size=1)) + 
  theme_minimal() +
  theme(legend.position = "none") +
  geom_vline(xintercept = 0, linetype = 'dotted') +
  geom_hline(yintercept = 0, linetype = 'dotted') +
  ggtitle(tf)

Version Author Date
182f9de GinoBonazza 2024-08-19
61e8d17 GinoBonazza 2024-08-01
# Run ulm
sample_acts <- run_ulm(mat=counts_mPAP_HC, net=net, .source='source', .target='target',
                  .mor='mor', minsize = 5)

# Transform to wide matrix
sample_acts_mat <- sample_acts %>%
  pivot_wider(id_cols = 'condition', names_from = 'source',
              values_from = 'score') %>%
  column_to_rownames('condition') %>%
  as.matrix()

# Get top tfs with more variable means across clusters
tfs <- arrange(f_contrast_acts, desc(score))$source
sample_acts_mat <- sample_acts_mat[,tfs]

# Scale per sample
sample_acts_mat <- scale(sample_acts_mat)
sample_acts_mat <- t(sample_acts_mat)
sample_acts_mat <- sample_acts_mat[, metadata_mPAP_HC$Sample, drop = FALSE]


# Choose color palette
palette_length = 100
my_color = colorRampPalette(c("Darkblue", "white","red"))(palette_length)

my_breaks <- c(seq(-3, 0, length.out=ceiling(palette_length/2) + 1),
               seq(0.05, 3, length.out=floor(palette_length/2)))

# Plot
pheatmap::pheatmap(sample_acts_mat, 
         border_color = NA, 
         color=my_color, 
         breaks = my_breaks,
         cluster_rows = F,
         cluster_cols = F,
         annotation_col = dplyr::select(metadata_mPAP_HC, -Batch_processing)) 

Version Author Date
61e8d17 GinoBonazza 2024-08-01
saveRDS(CM, 
        here::here(output_dir_data, "CM.rds"))

sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
 [3] LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8       
 [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8      
 [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
 [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8     
[11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] grid      stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] htmltools_0.5.8.1           DT_0.32                    
 [3] ComplexHeatmap_2.18.0       UpSetR_1.4.0               
 [5] ezRun_3.19.1                Biostrings_2.70.2          
 [7] XVector_0.42.0              enrichplot_1.22.0          
 [9] dorothea_1.14.1             OmnipathR_3.10.1           
[11] decoupleR_2.9.7             EnhancedVolcano_1.20.0     
[13] ggrepel_0.9.5               speckle_1.2.0              
[15] edgeR_4.0.16                limma_3.58.1               
[17] ReactomePA_1.46.0           AnnotationHub_3.10.0       
[19] BiocFileCache_2.10.1        dbplyr_2.4.0               
[21] org.Hs.eg.db_3.18.0         AnnotationDbi_1.64.1       
[23] clusterProfiler_4.10.1      gprofiler2_0.2.3           
[25] harmony_1.2.0               clustree_0.5.1             
[27] ggraph_2.2.1                rhdf5_2.46.1               
[29] Rcpp_1.0.12                 data.table_1.15.2          
[31] plyr_1.8.9                  gtable_0.3.4               
[33] gtools_3.9.5                ArchR_1.0.2                
[35] statmod_1.5.0               patchwork_1.2.0            
[37] scater_1.30.1               scuttle_1.12.0             
[39] Matrix.utils_0.9.7          RColorBrewer_1.1-3         
[41] scales_1.3.0                knitr_1.45                 
[43] gridExtra_2.3               png_0.1-8                  
[45] pheatmap_1.0.12             reshape2_1.4.4             
[47] lubridate_1.9.3             forcats_1.0.0              
[49] stringr_1.5.1               purrr_1.0.2                
[51] tidyr_1.3.1                 tibble_3.2.1               
[53] tidyverse_2.0.0             magrittr_2.0.3             
[55] ggplot2_3.5.0               dplyr_1.1.4                
[57] scCustomize_2.1.2           scDblFinder_1.16.0         
[59] Matrix_1.6-5                DropletUtils_1.22.0        
[61] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
[63] Biobase_2.62.0              GenomicRanges_1.54.1       
[65] GenomeInfoDb_1.38.7         IRanges_2.36.0             
[67] S4Vectors_0.40.2            BiocGenerics_0.48.1        
[69] MatrixGenerics_1.14.0       matrixStats_1.2.0          
[71] SeuratObject_5.0.2          Seurat_4.4.0               
[73] xlsx_0.6.5                  readxl_1.4.3               
[75] readr_2.1.5                 here_1.0.1                 

loaded via a namespace (and not attached):
  [1] igraph_2.0.3                  graph_1.80.0                 
  [3] ica_1.0-3                     plotly_4.10.4                
  [5] rematch2_2.1.2                zlibbioc_1.48.0              
  [7] tidyselect_1.2.1              rvest_1.0.4                  
  [9] bit_4.0.5                     doParallel_1.0.17            
 [11] clue_0.3-65                   lattice_0.22-5               
 [13] rjson_0.2.21                  blob_1.2.4                   
 [15] S4Arrays_1.2.1                parallel_4.3.1               
 [17] cli_3.6.2                     ggplotify_0.1.2              
 [19] goftest_1.2-3                 BiocIO_1.12.0                
 [21] bluster_1.12.0                grr_0.9.5                    
 [23] BiocNeighbors_1.20.2          uwot_0.1.16                  
 [25] shadowtext_0.1.3              curl_5.2.1                   
 [27] mime_0.12                     evaluate_0.23                
 [29] tidytree_0.4.6                leiden_0.4.3.1               
 [31] stringi_1.8.3                 backports_1.4.1              
 [33] XML_3.99-0.16.1               httpuv_1.6.14                
 [35] paletteer_1.6.0               rappdirs_0.3.3               
 [37] splines_4.3.1                 logger_0.3.0                 
 [39] bcellViper_1.38.0             sctransform_0.4.1            
 [41] ggbeeswarm_0.7.2              DBI_1.2.2                    
 [43] HDF5Array_1.30.1              jquerylib_0.1.4              
 [45] reactome.db_1.86.2            withr_3.0.0                  
 [47] git2r_0.33.0                  rprojroot_2.0.4              
 [49] xgboost_1.7.7.1               lmtest_0.9-40                
 [51] ggnewscale_0.4.10             tidygraph_1.3.1              
 [53] rtracklayer_1.62.0            BiocManager_1.30.22          
 [55] htmlwidgets_1.6.4             fs_1.6.3                     
 [57] labeling_0.4.3                SparseArray_1.2.4            
 [59] cellranger_1.1.0              reticulate_1.35.0            
 [61] zoo_1.8-12                    timechange_0.3.0             
 [63] foreach_1.5.2                 fansi_1.0.6                  
 [65] ggtree_3.10.1                 R.oo_1.26.0                  
 [67] irlba_2.3.5.1                 ggrastr_1.0.2                
 [69] gridGraphics_0.5-1            ellipsis_0.3.2               
 [71] lazyeval_0.2.2                yaml_2.3.8                   
 [73] survival_3.5-8                scattermore_1.2              
 [75] BiocVersion_3.18.1            crayon_1.5.2                 
 [77] RcppAnnoy_0.0.22              progressr_0.14.0             
 [79] tweenr_2.0.3                  later_1.3.2                  
 [81] ggridges_0.5.6                codetools_0.2-19             
 [83] GlobalOptions_0.1.2           KEGGREST_1.42.0              
 [85] Rtsne_0.17                    shape_1.4.6.1                
 [87] Rsamtools_2.18.0              filelock_1.0.3               
 [89] pkgconfig_2.0.3               xml2_1.3.6                   
 [91] GenomicAlignments_1.38.2      aplot_0.2.2                  
 [93] spatstat.sparse_3.0-3         ape_5.7-1                    
 [95] viridisLite_0.4.2             xtable_1.8-4                 
 [97] highr_0.10                    httr_1.4.7                   
 [99] tools_4.3.1                   globals_0.16.3               
[101] beeswarm_0.4.0                checkmate_2.3.1              
[103] nlme_3.1-164                  selectr_0.4-2                
[105] HDO.db_0.99.1                 crosstalk_1.2.1              
[107] digest_0.6.35                 farver_2.1.1                 
[109] tzdb_0.4.0                    yulab.utils_0.1.4            
[111] viridis_0.6.5                 glue_1.7.0                   
[113] cachem_1.0.8                  polyclip_1.10-6              
[115] generics_0.1.3                parallelly_1.37.1            
[117] ScaledMatrix_1.10.0           pbapply_1.7-2                
[119] vroom_1.6.5                   spam_2.10-0                  
[121] gson_0.1.0                    dqrng_0.3.2                  
[123] utf8_1.2.4                    graphlayouts_1.1.1           
[125] shiny_1.8.0                   GenomeInfoDbData_1.2.11      
[127] R.utils_2.12.3                rhdf5filters_1.14.1          
[129] RCurl_1.98-1.14               memoise_2.0.1                
[131] rmarkdown_2.26                R.methodsS3_1.8.2            
[133] future_1.33.1                 RANN_2.6.1                   
[135] Cairo_1.6-2                   spatstat.data_3.0-4          
[137] rstudioapi_0.15.0             cluster_2.1.6                
[139] whisker_0.4.1                 janitor_2.2.0                
[141] spatstat.utils_3.0-4          hms_1.1.3                    
[143] fitdistrplus_1.1-11           munsell_0.5.0                
[145] cowplot_1.1.3                 colorspace_2.1-0             
[147] rlang_1.1.4                   DelayedMatrixStats_1.24.0    
[149] sparseMatrixStats_1.14.0      dotCall64_1.1-1              
[151] ggforce_0.4.2                 circlize_0.4.16              
[153] xfun_0.42                     iterators_1.0.14             
[155] abind_1.4-5                   GOSemSim_2.28.1              
[157] interactiveDisplayBase_1.40.0 treeio_1.26.0                
[159] rJava_1.0-11                  Rhdf5lib_1.24.2              
[161] bitops_1.0-7                  promises_1.2.1               
[163] scatterpie_0.2.1              RSQLite_2.3.5                
[165] qvalue_2.34.0                 fgsea_1.28.0                 
[167] DelayedArray_0.28.0           GO.db_3.18.0                 
[169] compiler_4.3.1                prettyunits_1.2.0            
[171] beachmat_2.18.1               graphite_1.48.0              
[173] listenv_0.9.1                 workflowr_1.7.1              
[175] BiocSingular_1.18.0           tensor_1.5                   
[177] MASS_7.3-60.0.1               progress_1.2.3               
[179] BiocParallel_1.36.0           spatstat.random_3.2-3        
[181] R6_2.5.1                      fastmap_1.1.1                
[183] fastmatch_1.1-4               vipor_0.4.7                  
[185] ROCR_1.0-11                   rsvd_1.0.5                   
[187] KernSmooth_2.23-22            miniUI_0.1.1.1               
[189] deldir_2.0-4                  bit64_4.0.5                  
[191] spatstat.explore_3.2-6        lifecycle_1.0.4              
[193] ggprism_1.0.4                 restfulr_0.0.15              
[195] xlsxjars_0.6.1                sass_0.4.9                   
[197] vctrs_0.6.5                   spatstat.geom_3.2-9          
[199] snakecase_0.11.1              DOSE_3.28.2                  
[201] scran_1.30.2                  ggfun_0.1.4                  
[203] sp_2.1-3                      future.apply_1.11.1          
[205] bslib_0.6.1                   pillar_1.9.0                 
[207] metapod_1.10.1                locfit_1.5-9.9               
[209] jsonlite_1.8.8                GetoptLong_1.0.5